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Abstract 

Electrical transmission signals have been used for decades to characterize the internal structure of composite materials. We 
theoretically analyze the transmission of an electrical signal through a composite material which consists of two phases with 
different chemical compositions. We assume that the temperature of the biphasic system increases as a result of Joule 
heating and its electrical resistivity varies linearly with temperature; this last consideration leads to simultaneously study the 
electrical and thermal effects. We propose a nonlinear conjugate thermo-electric model, which is solved numerically to 
obtain the current density and temperature profiles for each phase. We study the effect of frequency, resistivities and 
thermal conductivities on the current density and temperature. We validate the prediction of the model with comparisons 
with experimental data obtained from rock characterization tests. 
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Introduction 

A composite medium can be defined as that made of at least two 
phases of different chemical compositions [1]. The study of 
composite media are of great interest to various areas such as 
physics, chemistry and materials science, among others [1-3]. 

The response of composite media when transmitting-absorbing 
waves of different intensities and frequencies has been analyzed by 
many previous studies. In [4], a fiber reinforced epoxy matrix 
composite is studied as electromagnetic wave absorbing material in 
a wide frequency range. Carbon fiber reinforced concrete [5] and 
metallic wire structures [6] have been characterized in terms of 
their capacity to absorb electromagnetic fields. It is interesting to 
note that some materials have the ability to shield electromagnetic 
waves, examples of these are styrene-butadiene rubber composites 
[7] , wood-cement boards [8] and nanocomposites [9] . 

The prediction of the effective properties of a composite 
medium is commonly based on the knowledge of both the volume 
fraction distribution and the value of the property of each of its 
phases [1] . A series of mathematical models have been developed 
to define global properties of a composite medium. In [1,10] the 
effective electrical conductivity is obtained from the individual 
electrical properties, also, heat transfer studies at macroscopic and 
microscopic levels have been conducted to determine the effective 
thermal conductivity [11,12]. To our knowledge the interaction 
between the heat diffusion and transmission of electric current 
through composite media has not been studied fuUy to date. 

The present theoretical model considers two non-deformable 
phases which are modeled as continua, and it is based on a 
previous model developed by Chavez and Mendez [13] who 
analyzed the conjugate heat and electromagnetic transfer mech- 



anism in a bimetallic conductor. Our system consists of a cylindral 
external phase confining a second cylindrical internal phase, 
henceforth external and internal phases, as depicted schematically 
in Figure 1. 

Maxwell's equations are coupled with the heat conduction 
equation considering a heat source term to account for Joule's 
effect. This coupling is performed for each phase; also, both phases 
are coupled with each other through the boundary conditions at 
the common interface. 

Particularly, the proposed model is intended to explain the role 
of some effects that occur during the electrical conduction 
processes in composite media, which to our knowledge have not 
been fuUy addressed: 

1 . The so-called skin effect observed in power transmission lines 
[14]. 

2. The Joule heating effect [15]. 

3. The effect of frequency on the current density and temperature 
distribution. 

4. The effect of volumetric fraction (porosity) on the bulk 
electrical resistivity. 

Modelling 

The physical model under study is a composite medium like that 
shown in Figure 1. We consider an external phase with radius b 
and an internal one with radius a. A sudden alternating electric 
current through this biphasic conductor is established. Thus, a rise 
of temperature results from the flow of electric current, caused by 
Joule's effect. We assumed that the resistivity in both phases varies 
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Figure 1. Schematic representation of a composite medium. 

doi:10.1371/journal.pone.0097895.g001 



linearly with temperature. For large values of the frequency 
associated to the alternating current, a redistribution of the current 
density is inevitable and the skin effect yields a tendency of the 
electric current to flow over the outer surfaces of both phases. In 
real conditions, the current density distribution depends on the 
values of electrical resistivities of each phase and the heat 
generated by Joule's effect which is transferred by heat conduction. 
Finally, we model the heat transferred to the environment by a 
convection process. 

Electromagnetic model 

Using Maxwell's equations, we can readily derive a wave 
equation to analyze the electromagnetic propagation. Therefore, 
the current density is governed by the following equation [16]: 



8J 

+y- 

dt 



(1) 



d^J, / 2(1) 8T l\ dJ, 
~dr^ \ \+(i>{T-T^,)~dr ^~r)~dr 
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where co is the frequency of the electrical signal, r is the radial 
coordinate, T is the temperature, is the environment 

temperature, /, is the current density function depending only 
on radial coordinate, ^ is the temperature coefficient for resistivity 
and / = \/— 1. 

In most practical cases, the term ycu^/x is smaller than iCO^i/A 
and can be neglected as a first approximation. In addition, we also 
introduce the well-known conductor skin depth parameter, 5, 
defined by (5 = (21/^ico)'^^ [17]. Thus, Eq. (2) can be rewritten for 
the internal phase as: 
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and for the external phase. 
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where the subscript "i" is used to denote the spatial dependence 
and the subscripts and "/" are used to denote external and 
internal phases, respectively. The "oo" subscript refers to external 
environmental conditions. 

The above equations system must be solved considering the 
following boundary conditions: 



where, I is the electric resistivity, / is the current density, fi is the 
magnetic permeability, y is the electric permitivity and t is the 
physical time. 

We consider only variations of the current density in the radial 
direction and the alternating current behaves like a sinusoidal 
wave. Therefore, the current density can be written as 
J = Js(r)e'"". On the other hand, the electrical resistivity has a 
linear variation with temperature [16] which can be written as 
/l = /loo [1 +<^(r— 7^00 )], and introducing it into Eq. (1) we obtain 
that, 



at r = 0 



dJsj 
dr 



atr = a: Aa),/[1 +<^/(7'/- 7"=o)l-/.!,/ = 
^-^,eI^+<I>e(Te-Too)]J,,e, 



(5) 



(6) 



(7) 
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at the surface r = b . 



-Jb- 



(8) 



f = 0 : 



(17) 



The boundary condition at the centre (Eq. (5)) is the symmetry 
condition, while the continuity of the electric field at the interface 

is expressed by Eq. (6). Eq. (7) refers to the continuity of the 
magnetic field, while Eq. (8) expresses a characteristic current 
density. Jh is the current density at the outer surface of the external 
phase and should be determined with the following restriction: 



I = 2n 



Jsjrdr + 



Js,Erdr 



(9) 



Thermal model 

The general heat diffusion equation can be expressed as [18]: 



In the above equations, k is the thermal conductivity, p is the 
density, c is the specific heat, h is the convective heat transfer 
coeHicient. 

Dimensional Analysis 

In order to reduce the number of physical parameters, we can 
perform a dimensional analysis. We first identify the characteristic 
convective time scale tc = (pc)ib/h. On the other hand, the 
suitable spatial scale is chosen as the radius of the external phase, 
r = b. Furthermore, the characteristic temperature drop ATc can 
be obtained through an energy balance between the heat 
generation term and the transient term, i.e.: 



WikVT]+q=ipc) 



dT 



(10) 



With this equation we can determine the gradients of 
temperature in both phases by taking into account their thermal 
properties and the amount of heat generated on each of them. 

Equation (10) is simplified by considering only temperature 
variations in the radial direction and exclusively heat generation 
by Joule's effect. Therefore, for the internal phase we have: 



AT, 



kpBi 



(18) 



where J a is the current density at the surface of the internal phase 
and Bi is the Biot number which measures the environmental 
conditions and is defined as 



Bi-- 



hb 
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and for the external phase: 



With the above set of characteristic geometrical and physical 
scales, the electromagnetic and thermal models can be simplified 
by introducing the following dimensionless variables and param- 
eters: 
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For the outer surface of the external phase, we have 



r = b: -kE^-^=h{TE-T^). (16) 



where <1> is the volumetric fraction defined as the ratio of the 
volume of the internal phase and the total volume of the composite 
media, while G is defined as the actual lenght of the inner path 
divided by the straight-line distance between the ends of the inner 
path, for porous media applications this parameter is known as 
tortuosity and is always larger or equal than one. 



Also an initial condition is necessary. Here we ha\'c considered 
that the composite media is initially at ambient temperature: 



Dimensionless electromagnetic model 

The system (3)-(8) can be rewritten in dimensionless form by 
using the above dimensionless parameters and variables. 
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The internal phase dimensionless model is: 
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and external phase dimensionless model is: 
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The corresponding dimensionless boundary conditions are: 
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Dimensionless thermal model 

In the same manner, we can use the dimensionless variables and 
parameters in order to obtain the following dimensionless thermal 
model: 

Therefore, the internal phase equation can be written as: 



\ d ( ^dd,\^kE{\+K,6,) 2 kEde, 



and external phase as: 



1 5 (jeE\ , XE(\+KEeE)^, ,2 



(pc)e ,.-SGe 



(26) 



(27) 



The associated boundary and initial conditions are given by: 

ddi 
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(30) 



(31) 



(32) 



The system of equations to be solved is formed by Eqs. (20) and 
(21) for the current density distribution coupled through the 
boundary conditions (23) and (24); and Eqs. (26) and (27) for the 
thermal behavior also coupled through their boundary conditions 
(29) and (30). It should be noticed that there is a coupling between 
the electromagnetic and thermal model due to the dependence of 
resistivity with temperature which is expressed by k parameter, 
henceforth the coupling parameter. 

Solution Methodology 

The above dimensionless electromagnetic and heat conduction 
equations, together with their boundary and initial conditions, 
represented here by the system of Eqs. (20)-(32) was solved by 
using a conventional iterative finite-differences method [19]. 

The electromagnetic equations are given by the system of Eqs. 
(20)-(25). These equations are complex because the right-hand 
sides of Eqs. (20) and (21) include, as a factor, the imaginary 
number /. Therefore, we separate for each region the electrical 
current density in a real part, </>^, and an imaginary part , 
through the relationship (p = <p^ + ig>'. The resulting equations are 
discretized together with the boundary conditions (22)-(25) 
considering central differences. In this form, we can construct a 
matrix system which can be solved with a simple Gauss 
elimination method. 

The corresponding equations for the thermal model given by 
the set of Eqs. (26) and (27) together with the boundary and initial 
conditions (28)-(32), require a different treatment. In this case, the 
above equations represent a non-stationary problem. Therefore, 
the numerical procedure is based on the well-known Crank- 
Nicholson finite difference scheme. In this manner, we obtain a 
tridiagonal matrix which is solved by the tridiagonal matrix 
algorithm (TDMA), also known as the Thomas algorithm. 

Finally, we introduce the following iterative scheme: firstly, a 
uniform profile for the temperature is considered, then we solve for 
the electrical ( urrcut dtnisity. In this manner, we can obtain the 
modulus or absolute value of this function. Introducing the above 
result into Eqs. (26) and (27), we obtain the first nonuniform 
temperature profile. Again, we can obtain a new current density 
and the foregoing procedure is repeated until a convergence 
criterion is fulfilled. This criterion is based on the comparison of 
the temperature and current density profiles. 

Validation 

It is important to note that if the coupling parameter is equal to 
zero (k = 0), the equations associated with the electromagnetic 
behavior, Eqs. (20) and (21), are no longer affected by the 
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temperature; thus, the system can be solved analytically. The 
solutions for the current density distributions, fiic) and (p,.(0 
(liquid and rock phases, respectively) are; 



and 



(33) 



(34) 



where Jo,J\, Yo, Yi denote the Bessel functions of first and second 
kind and of zeroth and first order, respectively. The variables 
/,A,Y,n,E and 0 are defined as: 



(l-0%/o l-i , l-i 

f= 7^' S= , h= 7==^, 35) 



This analytical solution has great relevance since it serves as a 
validation test for our numerical results. As shown in Fig. 2, the 
numerical simulations (open symbols) agree very well with the 
analytical solutions (solid lines). 

Results and Discussion 

In the system of equations (26)-(32) we have two couphng 
parameters Kj and K£. They, however, depend on each other 
because they are affected by the same ATc, thus 



ke= Ki- 



rn 



therefore it is necessary to know only one parameter. We choose 
the level of coupling between thermal and electrical model as 

K=Kj. 

In the same manner the skin parameter for the internal phase, 
Si, is related to that of the external one, because both phases 
transmit a wave at the same frequency. Thus we have: 



A = Ji{f)Yo(-f) + Jo{f)Yi(-f), 



(36) 



£E = SI- 



l-(0/G)2 V ^^''^^ 



(42) 



r = ,,iiXV^-l)(Jo(fVi(g)Yo(-h)-Mh)Ji(g)Yo(-f)), (37) 



n = Ji(f)Yoi-h) + Jo{h)Yi(-f), 



(38) 



E = (V^- l)8,fi,Ja(f)Ji (g) + V^EitiiMg)Ji if), (39) 



0 = ^siiJiMg) Ti ( -/) - (To - \)ErliMg) Yo(f), (40) 



From now on, £/ wiU be simply called £, the depth of 
penetration of the electrical signal. Therefore, the main param- 
eters in this problem are: IeP-i, kE/ki, {pc)^/(pc)i, ^Ie/ ^ii,S; Bi, 
K and fl). 

To study the transmission of the electrical signal in a composite 
medium, we perform a parametric analysis based on the following 
parameters: E, which is a function of the frequency, the ratio of 
resistivities, A^/l/, and the ratio of thermal conductivities, kE/ki. 
For simplicity, all these calculations assuming K=l, <1> = 0.5, 
(pc)^/(pc)i = \, fi^/fii = l and G=l; also, we consider Bi=\, 
which assures that the heat is efficiently transferred by convection 
from the external phase to the environment. To assure steady state 
solutions we performed all calculations using T = 20. 




Figure 2. Dimensionless current density distribution f as a 
function of the radial coordinate ^ for both conducting phases 
and three different values of the skin parameter c. Solid lines are 
computed from the analytical solution, the open symbols represent 
numerical simulation results. 
doi:1 0.1 371/journal.pone.0097895.g002 



Effect of e 

To analyze the effect of the skin parameter on the current 
density and temperature, numerical results were obtained from 
three different values of £ (0.1, 0.5 and 10) and are presented in 
Fig. 3. The series of three temperature profiles shown in Fig. 3a 
represents the steady-state solution for each corresponding value of 
£. For lower values of the skin parameter (which correspond to 
electric signals with high frequency), the temperature profiles 
become more uniform in both phases; on the contrary, when 
grows a slightly parabolic temperature profile is exhibited in the 
internal phase and even a steeper parabolic profde is observed for 
the external one. This parabolic behavior comes from the Joule 
effect (source therms in equations (26-27)). Figure 3b shows the 
current density distribution as a function of the dimensionless 
radial coordinate for the same conditions shown in Figure 3a. 
Clearly, for small values of £, the skin effect becomes noticeable in 
both phases. In particular the current density in the internal phase 
is higher in regions close to the interface [z^ajb). On the opposite 
side, high values of £ show nearly constant current distribution in 
both phases, hence the electric power tends to be transmitted as 
direct current. 
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Figure 3. Spatial profiles of temperature and current density, (a) Dimensionless temperature distribution, and (b) dimensionless current 
density distribution as a function of the dimensionless radial coordinate for three different values of the skin parameter £. For XeI)-i = W and 
kE/k, = 0.\. 

doi:1 0.1 371 /journal.pone.0097895.g003 



Effect of Xe/Ii 

Figure 4a shows the dimensionles.s temperature and Figure 4b 
current density distributions for three different values of the ratio 
of resistivities is/i/ (2, 3 and 10). For smaller values of the Ae/ 
ratio, it is expected that most of the current wiU flow through the 
internal phase which also increases the temperature due to the 
heat generated by Joule effect. In constrast, when the ratio aeI 
tends to increase, the current density profiles in the external phase 
are redistributed at the outer surface showing a more pronounced 
skin effect and smaller current densities gradients. The shape of 
the temperature profiles does not appear to change substantially. 
The temperature profile in the external phase becomes slightly 
parabolic in shape and its negative slope is less pronounced 
indicating smaller temperature gradients. 

Effect of kE/kj 

Figure 5a shows the dimensionless temperature distribution for 
three cUflferent values of the ratio of thermal conductivities kE / kj 



(0.1, 1.0 and 10). High kE/kj ratios imply an excellent heat 
transfer by diffusion in the external phase compared with the 
internal one, consequendy the heat generated at the internal phase 
faces more resistance to be transferred to the external one, which 
leads to a higher temperature gradient, while the external phase 
transfers the energy to the environment optimally. For the case 
when kE/ki = \ both media behave thermally as a single phase. 
On the other hand, when the ratio kE/k: is small the internal 
phase transfer the heat generated in a better way and as a result we 
obtain a flatter profile than in the external phase. Figure 5b shows 
the current density distribution under the influence of thermal 
conductivities (above mentioned). The electromagnetic equations 
(20) and (21) are affected by temperature gradients, as a result, the 
higher temperature gradient the more affected the current density 
profile; for the external phase where the temperature gradient is 
basically the same, the current density is almost not affected. 
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Figure 4. Spatial profiles of temperature and current density, (a) Dimensionless temperature distribution, and (b) dimensionless current 
density distribution as a function of the dimensionless radial coordinate for three different values of the ratio of resistivities i^//!/. For c = 0,5 and 

kE/k, = 0.\. 

doi:1 0.1 371 /journal.pone.0097895.g004 



Experimental Validation 

The measurement of the electrical conductivity of porous media 
has applications in many areas of science and technology such as 
geothermal reservoir engineering and soil science and is particu- 
larly important in the oil and gas industry to estimate the amount 
of petroleum in reservoirs. The main mechanism for electrical 
transport in porous rocks and soils is electrical conduction through 
the water filling the pore space. Furthermore, the pore-scale 
geometrical features of rocks and soils have a significant effect on 
their bulk electrical conductivity [20]. Conversely, the electrical 
conductivity can be used to infer the characteristics of pore space 
in rocks and soils [21]. 

The previous developed model can be used to study electro- 
magnetic transmission through rocks; which are a type of natural 
composite medium with two interpenetrating, percolating phases: 
the pore and solid networks [2] . It is worthy of note that rocks are 
also clasified as porous media [22], 



For the present application we suppose that the solid phase 
consists of a rock matrbc forming a solid cylinder (external phase). 
The fluid phase, typically a brine solution, represents a static 
incompressible fluid fdling a pore (internal phase) which is 
embedded into the sohd phase. From now in advance, all the 
variables with subindex E will be rename with the subindex r 
which stands for rock domain while for / we will use /, liquid 
domain. 

The relationship that links the electrical resistivity of a rock to its 
porosity was empirically determined by Archie [23]. He proposed 
that the formation factor depends on porosity obeying an inverse 
power law: 



i? = (D 



(43) 



where O is the bulk porosity and F is the formation factor, defined 
as the resistivity of a porous medium completely saturated with a 
conductive liquid divided by the hquid resistivity. The parameter 
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Figure 5. Spatial profiles of temperature and current density, (a) Dimensionless temperature distribution, and (b) dimensionless current 
density distribution as a function of the dimensionless radial coordinate for three different values of the ratio of thermal conductivity kE/ki. For 
;.£/;., = 10 and £ = 0.5. 
doi:1 0.1 371 /journal.pone.0097895.g005 



m is estimated by a fitting regression analysis of experimental data 
and is often called the cementation exponent; it implicidy provides 
information about the pore structure. 

Since its publication in 1942, Archie's law has been widely used 
in the petroleum industry to characterize reservoirs. This fact is 
easily explained by two aspects. First of all, the formula is simple 
and compact, which are undoubtedly attractive features. Secondly, 
this law describes to some extent the experimental trends of data 
from many diflFerent types of rock samples [24,25]. Nevertheless, it 
should be noted that the relation between F and O data acquired 
from reservoirs commonly show considerable scattering [24-26], 
which strongly suggests a non-trivial relationship between these 
variables. Therefore, Archie's law is simply a crude approxima- 
tion. 

In view of the previous definition for F and from the set of Eqs, 
(20)-(25), the following expression is constructed: 



<S/G 



(44) 



In the previous equation <1) is known as porosity, which 
previously was defined as volumetric fraction and now the 
tortuosity, G, is calculated according to the expression reported 
by Maciej [27] 



G=l-iS/«(f), 



(45) 



where /J is a parameter which relates the porosity with tortuosity. 

Typically in a porous medium applications F is shown as 
function of O in log-log plots, showing a negative slope behavior 
(power-law dependence) [23]. It is worth pointing out that for the 
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Figure 6. Comparison between the matliematical model 
(dasKied dotted line) and experimental data (crosses symbols) 
taken from [28]. All numerical results were computed with: 

A,//; = l X \0\ £=1 X 10'' and fi = 0.95. 
doi:1 0.1 371/journal.pone.0097895.g006 



Figure 7. Comparison between the mathematical model 
(dashed dotted line) and experimental data (crosses symbols) 
taken from [29]. All numerical results were computed with: 

;.,./i, = 625, £=1 X 10" and /? = 0.88. 
doi:10.1371/journal.pone.0097895.g007 



present model we fix the values for the thermal and electrical 
properties of rock and brine, and also an initial value for C> (which 
is increased at a constant A<t>); now we calculate F from equation 
(44) for each value of <I>. 

Figures 6 and 7 show the formation factor as a function of 
porosity. The prediction of our model is compared directly with 
experimental data [28] and [29]. The numerical model repre- 
sented by the dashed dotted line is solved with typical values of 
electrical resistivity for both phases [30]. The bulk magnetic 
permeability of the rock phase is of the same order than that of the 
liquid [31]. Moreover, the skin parameter was taken E=l x lO'', 
which involves frequencies in the range from 1 to 1000 Hz and the 
penetration of the electromagnetic wave in the whole domain for 
which experiments are usually conducted [32]. A remarkable 
agreement is observed between the numerical results and the 
experimental data. The model accurately reproduces the typical 
decreasing dependence of F with <J> for the experimental data; 
clearly, both the experimental data and our model for the 
formation factor do not follow a power-law dependence on 
porosity as proposed by Archie [23]. 

Conclusions 

A nonlinear theoretical model that describes the combined 
effect of heat transfer and transmission of an electromagnetic wave 
through a composite medium was developed. The model reveals 
the influence of the properties of the composite medium and the 
electrical signal on the current density and temperature as a 
function of the radial coordinate. Basically three parameters were 



varied: the skin depth for which the well-known skin effect is 
clearly observed for electric signals of high frequency; the ratio of 
electrical resistivities, showing diflFerent current distribution and 
their effect on the temperature through Joule's eflFect and finally 
the ratio of thermal conductivities which shows that the current 
density is clearly affected by temperature gradients. 

Additionally, our model allows the calculation of the formation 
factor without the use of empirical relations, such as Archie's law. 
Now, the formation factor can be calculated directly if the physical 
properties of the porous media and the electric signal are known. 
A curve of the formation factor as a function of the porosity shows 
a similar trend than that found in experimental data for a rock 
porous medium, clearly the behavior is non-linear and is well 
predicted by our model in strong contrast with the linear log-log 
Archie's law. To our knowledge, such a theoretical model does not 
exist in the specialized literature. 
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